#ifndef Algebra_sgArpack_h
#define Algebra_sgArpack_h

#include "AlgebraModule.h"
#include "EigenSolver.h"

SG_NAMESPACE_OPEN

namespace Algebra
{
    template <typename TDataType>
    struct TArpackPrivate;

    namespace arpack
    {
        /**
         * BMAT specifies the type of the matrix B that defines the
         * semi-inner product for the operator OP.
         *  BMAT = 'I' -> standard eigenvalue problem A*x = lambda*x
         *  BMAT = 'G' -> generalized eigenvalue problem A*x = lambda*M*x
         */
        enum class BMat : int
        {
            /// B = 'I' -> standard eigenvalue problem A*x = lambda*x
            Identity,
            /// B = 'G' -> generalized eigenvalue problem A*x = lambda*B*x
            Generalized
        };

        /**
        c  WHICH   Character*2.  (INPUT)
        c          'LM' -> want the NEV eigenvalues of largest magnitude.
        c          'SM' -> want the NEV eigenvalues of smallest magnitude.
        c          'LR' -> want the NEV eigenvalues of largest real part.
        c          'SR' -> want the NEV eigenvalues of smallest real part.
        c          'LI' -> want the NEV eigenvalues of largest imaginary part.
        c          'SI' -> want the NEV eigenvalues of smallest imaginary part.
         */
        enum class Which : int
        {
            /// 'LA' - compute the NEV largest (algebraic) eigenvalues
            LargestAlgebraic,
            /// 'SA' - compute the NEV smallest (algebraic) eigenvalues.
            SmallestAlgebraic,
            /// 'LM' - compute the NEV largest (in magnitude) eigenvalues.
            LargestMagnitude,
            /// 'SM' - compute the NEV smallest (in magnitude) eigenvalues.
            SmallestMagnitude,
            /// 'LR' - compute the NEV largest (real part) eigenvalues.
            LargestReal,
            /// 'SR' - compute the NEV smallest (real part) eigenvalues.
            SmallestReal,
            /// 'LI' - compute the NEV largest (imaginary part) eigenvalues.
            LargestImaginary,
            /// 'SI' - compute the NEV smallest (imaginary part) eigenvalues.
            SmallestImaginary,
            // /// 'BE' - compute NEV eigenvalues, half from each end of the
            // /// spectrum.  When NEV is odd, compute one more from the
            // /// high end than from the low end.
            // both_ends
        };

        /**
         *
         * IPARAM(7) = MODE
         * On INPUT determines what type of eigenproblem is being solved.
         * Must be 1,2,3,4,5; See under \Description of dsaupd  for the
         * five modes available.
         * @see https://github.com/opencollab/arpack-ng/blob/master/SRC/dsaupd.f
         *
         * IPARAM(7) = MODE
         * On INPUT determines what type of eigenproblem is being solved.
         * Must be 1,2,3; See under \Description of znaupd  for the four modes available.
         * @see https://github.com/opencollab/arpack-ng/blob/master/SRC/znaupd.f
         *
         */

        enum class Mode : int
        {
            MODE_1 = 1,  // mode 1
            MODE_2 = 2,  // mode 2
            MODE_3 = 3,   // mode 3
            MODE_4 = 4   // mode 4
        };
    }

    /**
     * @brief Wrap ARPACK-NG to solve eigenvalue problems.
     *
     * ARPACK-NG is a collection of Fortran77 subroutines designed to solve large
     * scale eigenvalue problems.
     *
     * @see https://github.com/opencollab/arpack-ng
     *
     * @author nene
     * @date September, 2024.
     */
    template <typename TDataType>
    class ALGEBRA_EXPORT TArpack : public TEigenSolver<TDataType>
    {
    public:
        ///@name Type Definitions
        ///@{
        typedef /*typename*/ TDataType                 DataType;
        typedef /*typename*/ Matrix_t<TDataType>       MatrixType;
        typedef /*typename*/ SparseMatrix_t<TDataType> SparseMatrixType;
        typedef /*typename*/ Vector_t<TDataType>       VectorType;
        ///@}

        TArpack ();
        ~TArpack ();        

        struct Option
        {
            arpack::BMat bmat = arpack::BMat::Generalized;  // the type of the matrix B
            // For double, mode 2/mode 3 supported
            // For complex, mode 3 supported.
            arpack::Mode mode  = arpack::Mode::MODE_3;  // default to aprack mode 3
            std::string which = "LM";  // Specify which of the Ritz values of OP to compute.
            int         nev = 1;    // Number of eigenvalues of OP to be computed. 0 < NEV < N.
            DataType   sigma = 0.0;  // shifts for mode 3
        };

        /// the default option to control the eigenvalue solver.
        Option getDefaultOption ();

        /// Number of eigenvalues of OP to be computed. 0 < NEV < N.
        void setOption (const Option& opt);

        /**
         * @brief solve the generalized eigenvalue problem A*x = lambda*B*x.
         *
         * Mode 3:  K*x = lambda*M*x, K symmetric, M symmetric positive semi-definite.
         *          ===> OP = (inv[K - sigma*M])*M  and  B = M.
         *          ===> Shift-and-Invert mode
         *
         * ARPACK's Mode 4 seems un-stable to check out the tageted eigenvalues. What's more,
         * the shift value 'sigma' must be given non-zero value, namely, sigma != 0.
         *
         * Thus, we mimic the Buckling mode using Mode 3. That's
         *
         * Mode 4:  K*x = lambda*KG*x, K symmetric positive semi-definite,
         *          KG symmetric indefinite
         *          ===> OP = (inv[KG - sigma*K])*K  and  B = K.
         *
         * @note the resultant eigenvalues will be presented in \f$ 1/lamda \f$ ascending algebraic order.
         *
         * @return [eigenvectors, eigenvalues]:
         *            eigenvectors = V(:. 0:nev),
         *            eigenvalues = V(0:nev, nev)
         */
        void solve (const SparseMatrixType&  A,
                    const SparseMatrixType&  B,
                    std::vector<DataType>&   eigenvalues,
                    std::vector<VectorType>& eigenvectors) override;

        // /// A matlab-like signature to solve the generalized eigenvalue problem
        // // see http://www.ece.northwestern.edu/local-apps/matlabhelp/techdoc/ref/eig.html
        // void solve (const SpMat_t& A, const SpMat_t& B, Mat_t& V, Mat_t& D);

        /// report whether previous computation was successful.
        Info_t info () const;

    private:
        // /// check the user options.
        // bool checkOptions ();

    private:
        TArpackPrivate<TDataType>* imp = nullptr;
    };

    // Enforce the explicit instantiation to hide the underlying implementation.

    // /**
    //  * dsaupd is usually called iteratively to solve one of the following problems:
    //  * 
    //  *   Mode 1:  A*x = lambda*x, A symmetric
    //  *            ===> OP = A  and  B = I.
    //  * 
    //  *   Mode 2:  A*x = lambda*M*x, A symmetric, M symmetric positive definite
    //  *            ===> OP = inv[M]*A  and  B = M.
    //  *            ===> (If M can be factored see remark 3 below)
    //  * 
    //  *   Mode 3:  K*x = lambda*M*x, K symmetric, M symmetric positive semi-definite
    //  *            ===> OP = (inv[K - sigma*M])*M  and  B = M.
    //  *            ===> Shift-and-Invert mode
    //  * 
    //  *   Mode 4:  K*x = lambda*KG*x, K symmetric positive semi-definite,
    //  *            KG symmetric indefinite
    //  *            ===> OP = (inv[K - sigma*KG])*K  and  B = K.
    //  *            ===> Buckling mode
    //  * 
    //  *   Mode 5:  A*x = lambda*M*x, A symmetric, M symmetric positive semi-definite
    //  *            ===> OP = inv[A - sigma*M]*[A + sigma*M]  and  B = M.
    //  *            ===> Cayley transformed mode
    //  * 
    //  * For TArpack<double>, The converged Ritz values are always returned in ascending
    //  * algebraic order.
    //  *
    //  * WHICH   Character*2.  (INPUT)
    //  *         Specify which of the Ritz values of OP to compute.
    //  *
    //  *         'LA' - compute the NEV largest (algebraic) eigenvalues.
    //  *         'SA' - compute the NEV smallest (algebraic) eigenvalues.
    //  *         'LM' - compute the NEV largest (in magnitude) eigenvalues.
    //  *         'SM' - compute the NEV smallest (in magnitude) eigenvalues.
    //  *         'BE' - compute NEV eigenvalues, half from each end of the
    //  *                spectrum.  When NEV is odd, compute one more from the
    //  *                high end than from the low end.
    //  *
    //  * @see https://github.com/opencollab/arpack-ng/blob/master/SRC/dsaupd.f
    //  */
    // template class ALGEBRA_EXPORT TArpack<Real_t>;

    // /**
    //  * znaupd  is usually called iteratively to solve one of the following problems:
    //  * 
    //  *   Mode 1:  A*x = lambda*x.
    //  *            ===> OP = A  and  B = I.
    //  * 
    //  *   Mode 2:  A*x = lambda*M*x, M hermitian positive definite
    //  *            ===> OP = inv[M]*A  and  B = M.
    //  *            ===> (If M can be factored see remark 3 below)
    //  * 
    //  *   Mode 3:  A*x = lambda*M*x, M hermitian semi-definite
    //  *            ===> OP =  inv[A - sigma*M]*M   and  B = M.
    //  *            ===> shift-and-invert mode
    //  *            If OP*x = amu*x, then lambda = sigma + 1/amu.
    //  * 
    //  * For TArpack<Complex_t>, it seems that the converged Ritz values are returned in
    //  * a way which is different from TArpack<double>.
    //  *
    //  * WHICH   Character*2.  (INPUT)
    //  *         'LM' -> want the NEV eigenvalues of largest magnitude.
    //  *         'SM' -> want the NEV eigenvalues of smallest magnitude.
    //  *         'LR' -> want the NEV eigenvalues of largest real part.
    //  *         'SR' -> want the NEV eigenvalues of smallest real part.
    //  *         'LI' -> want the NEV eigenvalues of largest imaginary part.
    //  *         'SI' -> want the NEV eigenvalues of smallest imaginary part.
    //  *
    //  * @see https://github.com/opencollab/arpack-ng/blob/master/SRC/znaupd.f
    //  */
    // template class ALGEBRA_EXPORT TArpack<Complex_t>;

    // // After that, export the routines which has been explicitly specialized.
    // // Otherwise, the compiler will complain that these routines could not be found.
    // template <>
    // ALGEBRA_EXPORT TArpack<double>::TArpack ();
    // template <>
    // ALGEBRA_EXPORT TArpack<double>::~TArpack ();
    // template <>
    // TArpack<double>::Option ALGEBRA_EXPORT TArpack<double>::getDefaultOption ();
    // template <>
    // void ALGEBRA_EXPORT TArpack<double>::setOption (const Option& opt);
    // template <>
    // Info_t ALGEBRA_EXPORT TArpack<double>::info () const;
    // template <>
    // Mat_t ALGEBRA_EXPORT TArpack<double>::solve (const SpMat_t& A, const SpMat_t& B);

    // template <>
    // ALGEBRA_EXPORT TArpack<Complex_t>::TArpack ();
    // template <>
    // ALGEBRA_EXPORT TArpack<Complex_t>::~TArpack ();
    // template <>
    // TArpack<Complex_t>::Option ALGEBRA_EXPORT TArpack<Complex_t>::getDefaultOption ();
    // template <>
    // void ALGEBRA_EXPORT TArpack<Complex_t>::setOption (const Option& opt);
    // template <>
    // Info_t ALGEBRA_EXPORT TArpack<Complex_t>::info () const;
    // template <>
    // ComplexMat_t ALGEBRA_EXPORT TArpack<Complex_t>::solve (const ComplexSpMat_t& A, const ComplexSpMat_t& B);

    // // alias for compatibilities with the external modules
    // typedef TArpack<Real_t> sgArpack;

}  // namespace Algebra

SG_NAMESPACE_CLOSE

#endif  // Algebra_sgArpack_h